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I  wish  to  deeply  acknowledge  the  Air  Force  Office  of  Scientific  Research  and 
particularly  Dr.  Nachman  for  supporting  our  work.  During  two  years  and  nine  months  of 
supported  research  we  performed  the  comprehensive  investigation  of  the  optical  polariton 
modes  bound  to  one-dimensional  arrays  of  dielectric  particles.  The  lowest  frequency 
resonant  modes  in  one-dimensional  arrays  of  dielectric  spherical  particles  were  studied. 
We  investigated  whether  they  can  exist,  calculated  their  lifetimes  (quality  factors), 
characterized  their  propagation  and  interference.  The  investigation  was  performed  using 
the  multisphere  Mie  scattering  formalism.  We  found  this  problem  interesting  practically 
because  these  systems  can  be  used  to  manipulate  optical  energy  in  the  sub-wavelength 
scale,  guide,  transfer,  filter,  store  and  covert  optical  energy.  High  quality  factor  modes 
can  also  be  used  in  microlasers. 

The  results  of  our  research  are  summarized  in  seven  published  articles  [1-7]  and 
one  submitted  paper  [8]  available  through  the  Physics  Preprint  Archive.  They  were 
reported  as  invited  talks  in  International  Conferences.9"11 

According  to  the  guidance  criterion  bound  modes  can  exist  if  the  half  of  a 
resonant  mode  wavelength  exceeds  the  interparticle  distance.  If  this  criterion  is  satisfied 
the  optical  energy  cannot  leave  the  particle  chain  because  of  the  conflict  with  the  energy 
and  momentum  conservation  laws  for  photon  emission  (light  cone  constraints).  We 
investigated  this  criterion  numerically  and  demonstrated  that  propagating  modes  exist  if 
the  refractive  index  of  dielectric  particles  exceeds  2. 1-6 

To  investigate  the  efficiency  of  light  binding  to  the  particle  chains  of  various 
shapes  we  calculated  quality  factors  of  most  bound  mode  depending  on  the  array  shape 
(circular  2  or  linear  chains^"6).  We  demonstrated  both  analytically  and  numerically  that 
the  quality  factor  of  the  most  bound  mode  in  circular  array  depends  exponentially  on  the 
number  of  particles  in  the  chain,  while  in  linear  chain  this  quality  factor  increases  as  the 

third  power  of  the  number  of 
particles.  Disordering  crucially 
affects  the  mode  quality  factor  in 
circular  arrays  while  in  linear 
particle  chains  its  effect  is  less 
significant. 

We  demonstrated  that 
modes  possessing  the  highest 
quality  factor  represent 
interesting  example  of  so  called 
slow  light  modes.6  These  modes 
have  a  vanishing  group  velocity 
in  the  limit  of  the  infinite  length 
of  the  particle  chain  because  they 
are  located  in  the  top  edge  of  a 
guiding  band.  Guiding  modes 
form  energy  bands  in  periodic 
chains  of  particles.  The  modes  in  the  upper  band  edge  have  a  zero  group  velocity  so  they 
behave  like  slow  light  modes. 

In  addition  to  quality  factor  we  investigated  the  propagation  of  guiding  modes. 
Wc  demonstrated  that  if  the  mode  is  excited  by  the  point  source  located  near  the  edge  of 


Fig.  1 .  Traffic  circle  model  to  study  polariton  interference 
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the  particle  chain,  then  practically  all  emitted  energy  is  absorbed  by  this  mode  that  can 
propagate  along  the  chain  without  losses  if  the  source  frequency  belongs  to  the  guiding 
band. 

Finally  we  investigated  whether  guiding  modes  can  interfere  with  each  other 
similarly  to  electronic  waves.  We  investigated  interference  in  the  traffic  circle  arrays^ 
(see  Fig.  1)  and  demonstrated  that  it  behaves  similarly  to  the  interference  of  electronic 
currents  in  two  wires.  It  is  thus  possible  to  realize  the  Aharonov  Bohm  like  effects  in 
particle  waveguides. 

To  complete  all  goals  described  above  we  used  the  multisphere  Mie  scattering 
formalism,  which  permitted  us  to  efficiently  solve  frequency  domain  Maxwell  equations 
for  the  system  containing  the  large  number  of  dielectric  spheres.  All  calculations  were 
performed  using  our  own  codes  given  in  the  Appendix  section.  The  calculations  were 
performed  using  Scilab  programming  package.12 

The  research  has  been  performed  with  the  help  of  tw  o  Pi’s  graduate  students  Olga 
Samoylova  and  Gail  Blaustein  partially  supported  by  this  project.  Olga  Samoylova  is 
planning  to  defend  her  PhD  in  May  2009  based  on  her  achievements,  while  Gail 
Blaustein  is  currently  support  by  the  SMART  Fellowship  Program  and  works  on  the 
other  project  associated  with  DNA  optical  properties.  Some  work  was  made  with  the  help 
of  Pi’s  collaborators  Dr.  ITya  Polishchuk  (Max  Planck  Institute  for  Physics  of  Complex 
Systems,  Dresden,  Germany)  and  Michael  Gozman  (Russian  Research  Center 
“Kurchatov  Institute,”  Moscow,  RUssia)  under  partial  support  of  this  project  and  the 
Tulane  University  Research  and  Enhancement  Fund. 
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Appendix.  Codes  in  the  Scilab12  programming  language  used  for  our  calculations.  They  can 
be  easily  converted  to  Matlab. 

function  y=Diag_int_A(r,  thet,  fi,  m,  n,  mu,  nu) 

//This  code  calculates  vector  translation  coefficients  A  as  function  of  polar  coordinates  for  translation 
//  vector,  m,  n;  mu,  nu  -  projections  of  angular  momentum  for  new  and  old  spherical  vector  functions 
m=-m; 

//Calculation  of  prefactor 

xl  =  ((G)Am)*exp((%i)*(mu+m)*fi)*(2*nu+l)/(2*n*(n+l)); 

x2  =  gammaln(n-m+ 1  )+gammaln(nu-mu+ 1  )-gammaln(n+m+ 1  )-gammaln(nu+mu+ 1 ); 
xx  =  xl  *exp(x2); 

//disp(xx); 

//Calculation  of  sum 

qmax=min(n,  nu,  (n+nu-abs(m+mu))/2); 

//Getting  the  Final  result 
SummK); 

//disp(qmax); 
for  q=0:qmax 
p=n+nu-2*q; 

uu=ab_GW(m,  n,  mu,  nu,  p); 
yl=uu(l); 

y3=LegendreXu(p,  m+mu,  cos(thct)); 

//y2=legendre(p,  abs(m+mu),  cos(thet)); 

//cc=m+mu, 

//if  cc  <  0 

//y2-y2*((- 1  )Acc)*exp(gammaln(p+cc+ 1  )-gammaln(p-cc+ 1 )); 

//end 

y22=y3*besselh(p-t-l/2,  r)*yl  *sqrt(%pi/(2*r)); 

//disp(y2); 

Summ=Summ+y22;//*(- 1  )Acc; 
end 

y=xx*Summ; 

endfunction 


function  y=Diag_int  A_Four(k,  q,  a,  m,  n,  mu,  nu,  Nmax) 

//This  code  calculate  Fourier  transform  of  A  for  linear  chain  of  period  a 

Summ^O; 
for  ii=l  :Nmax 


FA9550-06-1  -01 1 0,  Optical  excitations  and  energy  transfer  in  nanoparticle  waveguides 
Final  performance  report 


A.  L.  Burin 


-4- 


Summ=Summ+exp(%i*k*a*ii)*Diag_int_A(a*ii*q,  0,  0,  m,  n,  mu,  nu)+exp(- 
%i*k*a*ii)*Diag_int_A(a*ii*q,  %pi,  0,  m,  n,  mu,  nu); 
end 

y=Summ; 

endfunetion 

function  y=Diag_int_A_Four_xy_sph(mm,  q,  a,  m,  n,  mu,  nu,  N) 

//This  eode  ealeulate  Fourier  transform  of  A  for  the  spherical  array  in  x-y  plane 

Summ=0; 

R~a*N/(2*%pi); 
for  ii=l  :N-1 

//  Summ=Summ+exp(%i*mm*ii*2*%pi/N)*Diagjnt_A(2*q*R*sin(ii*%pi/N),  %pi/2,  n*%pi/N,  m,  n, 
mu,  nu); 

Summ=Summ+exp(%i*mm*ii*2*%pi/N)*Diag_int_A(2*q*R*sin(ii*%pi/N),  %pi/2,  0,  m,  n,  mu,  nu); 
//  disp(Summ); 
end 

y=Summ; 

endfunetion 


function  y=ab_GW(m,  n,  mu,  nu,  p) 

//  Function  used  in  evaluation  of  diagonal  vector  translation  coefficients 
//  Calculation  of  diagonal  a-part 
//  and  off-diagonal  b-part 

z=(-l)A(m+mu)*(2*p+l)*Wigner3j(n,  nu,  p,  0,  0,  0)*Wigner3j(n,  nu,  p,  m,  mu,  -m-mu); 

zz=gammaln(n+m+ 1  )+gammaln(nu+mu+ 1  )+gammaln(p-m-mu+ 1  )-gammaln(n-m4  1  )-gammaln(nu-mu+ 1  )- 

gammaln(p+m+mu+ 1 ); 

zz=zz/2; 

yy=z*exp(zz); 

//Calculation  of  off-diagonal  b-part 

zl=(-l)A(m+mu)*(2*p+3)*Wigner3j(n,  nu,  p,  0,  0,  0)*Wigner3j(n,  nu,  p+1,  m,  mu,  -m-mu); 

zzl  =gammaln(n+m+ 1  )+gammaln(nu+mu+ 1  )+gammaln(p-m-mu+2)-gammaln(n-m+ 1  )-gammaln(nu- 

mu+ 1  )-gammaln(p+irri-mu+2); 

zzl=zzl/2; 

yyl=zl*exp(zzl); 

y(  1  )=yy*((%i)Ap)*(n  *(n+ 1 )+  nu*(nu+ 1  )-p*(p  *  1 )), 
y(2)=yy  1  *((%i)A(p+ 1  ))*sqrt(((p+ 1  )A2-(n-nu)A2)*((n+nu+ 1  )A2-(p+ 1  )A2)); 

endfunetion 

function  y=Wigner3j(j  1,  j2,  j3,  ml,  m2,  m3) 

//  evaluation  of  WIgner  3j-symbols  used  in  vector  translation  coefficients 

if  (ml+m2+m3  — 0) 

z=0; 

else 

kmax  =  min(j  Hj2-j3,  jl-ml,  j2+m2); 
kmin  -  max(0,  j2-j3-ml,  j  l-j3+m2) 
zl=(-l)A(jl+j2+m3); 

z2  =  gammaln(jl-ml+l)+gammaln01+ml  +  l)+gammaln(j2-m2+l)+gammaln(j2+m2-+ 1)+  gammaln(j3- 
m3+ 1  )+gammaln(j3+m3+ 1 ); 

z2=z2-gammaln(j  1  +j2-j3+ 1  )-gammaln(j  1  -j2+j3+ 1  )-gammaln(-j  1  +j2+j3+ 1  )-gammaln(j  1  +j2+j3+2); 
/2=z2/2; 
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Summ=0; 

for  k=kmin:kmax 

z3  =  gamma!  n(j  1  +j2-j3+ 1  )+gammaln(j  1  -j2+j3+  1  )+gammaln(-j  1  +j2+j3+ 1  )-gammaln(k+ 1  )-gammaln(j  l+j2- 
j3-k+l); 

z3=z3-gammaln(j  1  -m  1  -k+ 1  )-gammaln(-j2+j3+ml  +k+ 1 )  -  gammaln(j2+m2-k+ 1  )-gammaln(-j  1  +j3- 
m2+k+l); 

tt=((- 1  )Ak)*exp(z3+z2); 

Summ=Summ+tt; 

end 

z=Summ*zl; 

end 

y=z 

endfunction 


function  y=LegendreXu(l,  m,  x) 

//Redefinition  of  Legendre  polynomials  in  accordance  with  the  multisphcre  Mie  scattering  formalism 
if  abs(x)~l 
if  (m>=0) 

z=((-l  )Am)*legendre(l,  m,  x); 
else 
m  =-m; 

z=((-l)Am)*legendre(l,  m,  x)*((-l)Am)*exp(gammaln(l-m+l  )-gammaln(l+m+l )); 

end 

else 

if  m— 0 

z=0; 

else 

if  x=l 

z=l; 

else 

z=(-l)Al; 

end 

end 

end 

y=z; 

endfunction 

function  y=LegendreXuDer(l,  m,  x) 

//Evaluation  of  of  Legendre  polynomial  derivatives 
if  abs(x)~l 
if  1=0 
zl=0; 
else 

if  abs(m)<l 

zl=(l*x*LegendreXu(l,  m,  x)-(l+abs(m))*LegendreXu(l-l,  m,  x))/sqrt(  1 -xA2); 
else 

zl=(l*x*LegendreXu(l,  m,  x))/sqrt(l-xA2), 
end 
end 
else 

if  abs(m)~l 

zl=0; 

else 

zl  =l*(l+l)/2; 

zl  =  zl*LegendreXu(l,  m,  l/2)/LegendreXu(l,  abs(m),  l/2)*(x)A(l+l ); 
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end 

end 

y=zi; 

endfunction 

//Alternative  approaches  to  vector  translation  coefficients  A 

function  y=Diag_int_A  l(r,  thet,  fi ,  m,  n,  mu,  nu) 
z  1  =((- 1  )A(n+m))*(2*n+ 1  )/(2*n*(n+ 1  ))*exp(%i  *(mu-m)*fi); 
qmax=min(n,  nu,  (n+nu-abs(m-mu))/2); 

Summ=0; 
for  q=0:qmax 
p=n+nu-2*q; 

Summ=Sumrn+Gauntl(-m,  n,  mu,  nu,  n+riu-2*q)*((-l)Aq)*(n*(n+l)-t-nu*(nu+l)- 
p*(p+l))*sqrt(%pi/(2*r))*besselh(p+l/2,  r)*legendrc(p,  mu-m,  cos(thet)); 
end 

y=Summ*zl ; 
endfunction 

function  y=Diagjnt_A2(r,  thet,  fi,  m,  n,  mu,  nu) 

z  1  =((- 1  )A(n+m))*exp(%i  *(mu-m)*fi)*sqrt((2*n+ 1  )*(2*nu+ 1  )/(n*(n+ 1  )*nu*(nu+ 1 ))); 
z2  =  (gammaln(n+m+ 1  )+gamrnaln(nu-mu+ 1  )-gammaln(n-m + 1  )-gammaln(nu+mu+l  ))/2; 
zl=zl  *exp(z2); 

qmax=min(n,  nu,  (n+nu-abs(m-mu))/2); 

Summ=0; 
for  q=0:qmax 
p^n+nu-2*q; 

Summ=Summ+Gauntl  (-m,  n,  mu,  nu,  n+nu-2*q)*((%i)Ap)*(n*(n+l  )+nu*(nu+l)- 
p*(p+l  ))*sqrt(%pi/(2*r))*besselh(p+l/2,  r)*legendre(p,  mu-m,  cos(thet)); 
end 

y=Summ*zl; 

endfunction 

function  y=OffDiag_int_B(r,  thet,  fl,  m,  n,  mu,  nu) 

//This  code  calculates  vector  translation  coefficients  B  as  function  of  polar  coordinates  for  translation 
//  vector,  m,  n;  mu,  nu  -  projections  of  angular  momentum  for  new  and  old  spherical  vector  functions 
m=-m; 

//Calculation  of  prefactor 

x  1  =  ((- 1  )A(m+ 1  ))*cxp((%i)*(mu+m)*fi)*(2*nu+ 1  )/(2*n*(n+ 1 )); 

x2  =  gammaln(n-m+ 1  )+gammaln(nu-mu+ 1  )-gammaln(n+m+ 1  )-gammaln(nu+mu+ 1 ); 

xx  =  xl  *exp(x2); 

//disp(xx); 

//Calculation  of  sum 

qmax=min(n,  nu,  (n+nu+l-abs(m+mu))/2); 

//Getting  the  final  result 
Summ=0; 

//disp(qmax); 
for  q=l  :qmax 
p=n+nu-2*q; 

//disp(p); 

uu=ab_GW(m,  n,  mu,  nu,  p); 
yl=uu(2); 

//disp(yl); 

y3=LegcndreXu(p+l ,  m+mu,  cos(thet)); 

//y2=legendrc(p+l,  abs(m+mu),  cos(thet)); 

//cc=m+mu; 
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//if  cc  <  0 

//y2=y2*((- 1  )Acc)*exp(gammaln(p+ 1  +cc+ 1  )-gammaln(p+ 1  -cc+ 1 )); 

//disp(y2); 

//end 

y22=y3*besselh(p+3/2,  r)*y  1  *sqrt(%pi/(2*r)); 

//disp(y2); 

Summ=Summ+y22;//*(- 1  )Acc; 
end 

y=xx*Summ; 

endfunction 

function  y=OffDiagjnt_B_Four_xy_sph(mm,  q,  a,  m,  n,  mu,  nu,  N) 

//This  code  calculate  Fourier  transform  of  B  for  the  spherical  array  in  x-y  plane 

Summ-0; 

R=a*N/(2*%pi); 
for  ii=l  :N-I 

Summ=Summ+exp(%i*mm*ii*2*%pi/N)*OffDiagjnt_B(2*q*R*sin(ii*%pi/N),  %pi/2,  0,  m,  n,  mu, 
nu), 
end 

y=Summ; 

endfunction 

function  y=Miel(n,  z,  n_r) 

//Mie  scattering  coefficient  for  electric  type  scattering 

num=(rieatti  besselj(n,  z)*rieatti_besselj_der(n,  n_r*/)-n_r*rieatti_besselj_der(n,  z)*neatti_besselj(n, 
n  r*z)); 

den=(ricatti_bcsselh(n,  z)*ricatti  besselj_der(n,  n_r*z)-n_r*ricatti_bcssclh_der(n,  z)*ricatti_bcssclj(n, 

n_r*z)); 

a_n=num/den; 

y-a_n 

endfunction 

function  y=Mie2(n,  z,  n_r) 

//Mie  scattering  coefficient  for  magnetic  type  scattering 

num=(ricatti  besselj(n,  z)*ncatti  besselj  der(n,  n_r*z)*n  r-ricatti_bcsselj_dcr(n,  z)*rieatti_besselj(n, 
n_r*z)); 

den=(ricatti_besselh(n,  z)*ncatti_besselj_der(n,  n_r*z)*n_r-ncatti_bcssclh_dcr(n,  z)*ncatti  bcssclj(n, 

n_r*z)); 

b_n=num/den; 

y^b_n 

endfunction 

//Functions  below  use  the  generalized  Newton-Raphson  algorithm  to  calculate  quality  factor  for  various 
//modes 

function  y=SolveMie2DipSph(n_r,  n,  m,  N,  mm,  q_0) 

//  Solver  for  dipolar  approach  and  lowest  Mie  resonance  for  closely  packed  circular  array 
//  n_r  -refractive  index  2.7  -  TiO_2,  3.5  -  GaAs,  2  -  ZnO 
//  n=l  -  dipoles,  m=0  -  tl,  m=l  - 12,  m=-l  - 1 
//  N  -  number  of  spheres 

//  mm=0,  1 , .  .N  angular  momentum  of  mode  mm=N/2  -  most  interesting  ease 
//  q_0  -  initial  value  for  iteration  procedure 

//  Assuming  the  distance  between  sphere  is  2,  targets  are  the  quality  factor  and  the  decay  rate  for  //a=200nm 

h =0.0000000001; 

q=q_0; 
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if  m=0 

disp("l  st  transverse  mode  tlH); 

delt  -  1 ; 

k=0; 

while  abs(delt/q)  >  hA(3/2) 
k=k+ 1 ; 

zz  =  1/Mie2(l,  q  ,n_r)  +  Diag_int_A_Four_xy_sph(mm,  q,  2,  0,  n,  0,  n,  N); 

//  disp(l/Mie2(l ,  q,  n_r)); 

//  disp(mm); 

//  disp(n); 

//  disp(N); 

//  disp(Diag_int_A_Four_xy  sph(mm,  q,  2,  0,  n,  0,  n,  N)); 

//  disp(Diag_int_A_Four_xy_sph(50,  q,  2,  0,  1,0,  1,  100)); 

//  disp(q); 

//  disp(zz); 

zz_h=l/Mie2(l,  q+h  ,n_r)  +  Diag_int_A_Four_xy_sph(mm,  q+h,  2,  0,  n,  0,  n,  N); 
zz_der  =  (zz_h-zz)/h; 
delt=-zz/zz  der; 
q=q+delt; 
if  (k>20) 
delt=0; 
end 
end 

disp(k); 

end 

if  m=l 

disp("2nd  transverse  mode  t2M); 

delt  =  1 ; 

k=0; 

while  abs(delt/q)  >  hA(3/2) 
k=k+ 1 ; 

Av=  (Diag _ int  A  Four _xy_sph(mm+ 1 ,  q,  2,  1,  n,  1,  n,  N)+Diag_int  A_Four  _xy  sph(mm-l,  q,  2,  1,  n, 

1,  n,  N))/2; 

Diff=(Diag_int_A_Four_xy_sph(mm+ 1 ,  q,  2,  l,n,  1,  n,  N)-Diag_int_A_Four_xy_sph(mm- 1 ,  q,  2,  1,  n, 
1,  n,  N))/2; 

Repuls  =  Diag_int_A_Four_xy_sph(mm,  q,  2,  1,  n,  -l,n,  N)*Diag_int_A_Four_xy_sph(mm,  q,  2,  -1,  n, 

1, n,  N); 

zz  =  1/Mie2(l ,  q,  n  r)  +  A v+sqrt(DiffA2+ Repuls); 

//  zz  =  1/Mie2(l ,  q,  n_r)+Av  +  Diag_int_A_Four_xy_sph(mm,  q,  2,  1 ,  n,  1 ,  n)/2; 

Av  h=  (Diag_int_A_Four_xy_sph(mm+l,  q+h,  2,  1,  n,  1,  n,  N)+Diag_int_A_Four_xy_sph(mm-l ,  q+h, 

2,  l,n,  1, n, N))/2; 

Diff  h=(Diag_int_A_Four_xy_sph(mm+l ,  q+h,  2,  1,  n,  1,  n,  N)-Diag_int_A_Four_xy_sph(mm-l ,  q+h, 
2,  1,  n,  1,  n,  N))/2; 

Repuls_h  =  Diag_int_A_Four  xy_sph(mm,  q+h,  2,  1,  n,  -1,  n,  N)*Diag_int_A_Four_xy  sph(mm,  q+h, 
2,  -1 ,  n,  1 ,  n,  N); 

zz_h  =  1/Mie2(l,  q+h,  n  r)  +  Av_h+sqrt(Diff_hA2+Repuls_h); 

//  zz_h  =  1/Mie2(l,  q+h,  n_r)+Av_h  +  Diagjnt_A_Four_xy_sph(mm,  q+h,  2,  1,  n,  1,  n)/2; 

//  zz  =  1/Mie2(l ,  q  ,n_r)  +  Diag_int_A_Four_xy_sph(mm,  q,  2,  1,  n,  1,  n, 

N)+sqrt(Diag_int_A  Four_xy_sph(mm,  q,  2,  1,  n,  -1,  n,  N)*Diag_int_A_Four_xy_sph(mm,  q,  2,  -1,  n,  1 , 
n,  N)); 

//  zz_h=l/Mie2(l,  q+h  ,n_r)  +  Diag_int_A_Four_xy_sph(mm,  q+h,  2,  1,  n,  1,  n, 

N)+sqrt(Diag  int_A  Four_xy_sph(mm,  q+h,  2,  1,  n,  -1,  n,  N)*Diag  int_A_Four_xy_sph(mm,  q+h,  2,  -1, 
n,  1,  n,  N)); 

//  zz  =  1/Mie2(l ,  q  ,n_r)  +  Diag_int_A_Four_xy_sph(mm,  q,  2,  1,  n,  1,  n,  N)+ 

Diag_int_A_Four_xy  sph(mm,  q,  2,  1,  n,  -1,  n,  N)/2; 
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//  zz_h=l/Mie2(l ,  q+h  ,n_r)  +  Diag_int_A_Four_xy_sph(mm,  q+h,  2,  1 ,  n,  1 ,  n, 
N)+Diag_int_A_Four_xy_sph(mm,  q+h,  2,  1,  n,  -1,  n,  N)/2; 


zz_der  =  (zz_h-zz)/h; 
dclt=-zz/zz_dcr; 
q=q+delt; 
if  k>=20 
delt=0; 
end 

//  disp(q); 

//  disp(zz); 

//  disp(delt); 
end 

disp(k); 

end 


if  m=-l 

dispC'longitudinal  mode  I"); 

delt  =  1 ; 

k=0; 

while  abs(delt/q)  >  hA(3/2) 
k=k+l; 

Av=  (Diag_int_A_Four_xy_sph(rnrn+l,  q,  2,  1,  n,  1,  n,  N)+Diag_int_A_Four_xy_sph(mm-l ,  q,  2,  1,  n, 

1 ,  n,  N))/2; 

Diff=(Diag  int  A  Four_xy_sph(mm+l,  q,  2,  l,  n,  l,  n,  N)-Diag_int_A_Four_xy_sph(mm-l ,  q,  2,  1,  n, 

1 ,  n,  N))/2; 

Repuls  =  Diag_int_A_Four_xy_sph(rnrn,  q,  2,  1,  n,  -1,  n,  N)*Diagjnt_A_Four_xysph(rnrn,  q,  2,  -1,  n, 
U  n,  N); 

zz  -  1/Mie2(l,  q,  n  r)  +  Av-sqrt(DiffA2+Repuls); 

Av_h=  (Diag_int_A_Four_xy_sph(mm+l,  q  +  h,  2,  1,  n,  1,  n,  N)+  Diag_int_A_Four_xy_sph(mm-l,  q+h, 

2,  1 ,  n,  1 ,  n,  N))/2; 

Diff_h=(Diag  int  A  Four  xy  sph(mm+l,  q+h,  2,  l,n,  1 ,  n,  N)-Diag_int_A_Four_xy_sph(mm-l ,  q+h, 
2,  1 ,  n,  1 ,  n,  N))/2; 

Repuls  h  =  Diag_int_A_Four_xy_sph(mm,  q+h,  2,  1,  n,  -1,  n,  N)*Diag_int_A_Four_xy_sph(mrn,  q+h, 
2,  -1,  n,  1 ,  n,  N); 

zz_h  =  1/Mie2(l,  q+h,  n  r)  +  Av_h-sqrt(Diff_hA2+Repuls_h); 
zz_der  =  (zz_h-zz)/h; 
delt=-zz/zz_der; 
q=q+delt; 
if  k>20 
delt=0; 
end 
end 

disp(k); 

end 

if  imag(q)~0 
disp(-real(q)/(2*imag(q))); 
end 

y=q; 

endfunction 

function  y=SolveMie2abSph(n_r,  N,  q_0) 

//  Solver  for  dipolar  approach  +  off-diagonal  interaction  and  lowest  Mie  resonance  for  closely  packed 
//  circular  array  taking  into  account  the  off-diagonal  interaction 
//  n  r  -refractive  index  2.7  -  TiO_2,  3.5  -  GaAs,  2  -  ZnO 
//  n=l  -  dipoles,  m=0  -  tl,  m=l  - 12,  m=-l  - 1 
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//  N  -  number  of  spheres 

//  mm=0,  1,  ...N  angular  momentum  of  mode  mm=N/2  -  most  interesting  case 
//  q_0  -  initial  value  for  iteration  procedure 

//  Assuming  the  distance  between  sphere  is  2,  targets  are  the  quality  factor  and  the  decay  rate  for 

//  a=200nm 

h=0.0000000001; 

q=q_0; 

n=l; 

mm=N/2; 

disp("  1  st  transverse  mode  1 1 "); 

delt  =  1 ; 

k=0; 

while  abs(delt/q)  >  hA(3/2) 
k=k+ 1 ; 

denomin=l/Miel(l,  q,  n_r)+Diag_int_A_Four_xy_sph(mm-l,  q,  2,  1,  n,  1,  n,  N)- 
2*Diag_int_A_Four_xy_sph(mm,  q,  2,  -1,  n,  1,  n,  N); 

zz  =  1/Mie2(l ,  q  ,n_r)  +  Diag_int_A_Four_xy_sph(mm,  q,  2,  0,  n,  0,  n,  N)- 
4*OffDiag_int  B_Four_xy_sph(rnm-l/2,  q,  2,  0,  1,  1,  1,  N)A2/dcnomin  ; 

denominh=l/Miel(l,  q+h,  n  r)+Diagjnt_A_Four_xy_sph(mm-l,  q+h,  2,  l,n,  1,  n,  N)- 
2*Diag_int_A_Four_xy_sph(mm,  q+h,  2,  -1,  n,  1,  n,  N); 

zz_h=l/Mie2(l ,  q+h  ,n_r)  +  Diag_int_A_Four_xy_sph(mm,  q+h,  2,  0,  n,  0,  n,  N)- 
4*OffDiag_int_B_Four_xy_sph(mm-l/2,  q+h,  2,  0,  1,  1,  1,  N)A2/dcnominh; 
zz_der  =  (zz_h-zz)/h; 
delt=-zz/zz  der; 
q=q+delt; 
if  (k>20) 
delt=0; 
end 
end 

disp(k); 

disp(-real(q)/(2*imag(q))); 

y=q; 

endfunction 


function  y=SolveMie2StrangeSph(n_r,  N,  q_0) 

//  Solver  for  dipolar  approach  +  off-diagonal  interaction  and  lowest  Mie  resonance  for  closely  packed 
//circular  array  taking  into  account  the  off-diagonal  interaction 
//  n_r  -refractive  index  2.7  -  TiO_2,  3.5  -  GaAs,  2  -  ZnO 
//  n=l  -  dipoles,  m=0  -  tl ,  m=l  -  t2,  m=-l  - 1 
//  N  -  number  of  spheres 

//  mm=0,  1 ,  ...N  angular  momentum  of  mode  mm=N/2  -  most  interesting  case 
//  q_0  -  initial  value  for  iteration  procedure 

//  Assuming  the  distance  between  sphere  is  2,  targets  are  the  quality  factor  and  the  decay  rate  for  //a=200nm 

h=0.0000000001; 

q=q_0; 

n=l ; 

mm=N/2; 

disp(nlst  transverse  mode  tl"); 

delt  =  1 ; 

k=0; 

while  abs(delt/q)  >  hA(3/2) 
k=k+ 1 ; 

dcnomin=l/Miel(l,  q,  n_r)+Diag_int_A_Four_xy_sph(mm-l ,  q,  2,  1,  n,  1,  n,  N)- 
2*Diag_int_A_Four_xy_sph(mm,  q,  2,  -1,  n,  1,  n,  N); 
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zz  =  1  /Mie2(  1 ,  q  ,n_r)  +  Diag_int_A_Four_xy_sph(mm,  q,  2,  0,  n,  0,  n,  N)- 
4*OffDiag_int_B_Four_xy_sph(mm-l/2,  q,  2,  0,  1,  1,  1,  N)A2/denomin  ; 

// bl-a2  00 

denominl  =  l/Miel(2,  q,  n_r)+Diag_int_A_Four_xy_sph(mrn,  q,  2,  0,  2,  0,  2,  N); 
zz=zz-OfFDiag_int_B_Four_xy_sph(mm,  q,  2,  0,  1 , 0,  2,  N)*OffDiag_int_B_Four_xy_sph(mm,  q,  2,  0, 

2,  0,  1 ,  N)/denominl ; 

//bl-b3  00 

//  denomin2=l/Mie2(3,  q,  n_r)+Diag_int_A_Four_xy_sph(mrn,  q,  2,  0,  3,  0,  3,  N); 

//  zz=zz-Diag_int_A_Four_xy_sph(mm,  q,  2,0,  1 , 0,  3,  N)*Diag_int_A_Four_xy_sph(mm,  q,  2,  0,  3,  0, 

1,  N)/denomin2; 

//bl  -a2  02 

denomin3^1/Miel(2,  q,  n_r)+Diag_int_A_Four_xy_sph(mm,  q,  2,  2,  2,  2,  2,  N); 
zz=zz-OffI)iag_int_B_Four_xy_sph(rnrn,  q,  2,  0,  1, 2,  2,  N)*OffDiag_int_B_Four_xy_sph(mm,  q,  2,  2, 

2,  0,  1,  N)/denomin3; 

//bl-a2  0-2 

denomin4=l/Miel(2,  q,  n_r)+Diag_int_A_Four_xy_sph(rnrn,  q,  2,  -2,  2,  -2,  2,  N); 
zz=zz-OffDiag_int_B_Four_xy_sph(mm,  q,  2,0,  1 ,  -2,  2,  N)*OffDiag_int_B_Four_xy_sph(mm,  q,  2,  -2, 
2,  0,  1,  N)/denomin4; 

//bl  -b3  02 

//  denomin5=l/Mie2(3,  q,  n_r)+Diag_int_A_Four_xy_sph(mm,  q,  2,  2,  3,  2,  3,  N); 

//  zz=zz-DiagJnt_A_Four_xy_sph(mm,  q,  2,  0,  1, 2,  3,  N)*Diag_int_A_Four_xy_sph(rnrn,  q,  2,  2,  3,  0, 

1 ,  N)/denomin5; 

//bl  -b3  0-2 

//  denomin6=l/Mie2(3,  q,  n_r)+Diag_int_A_Four_xy_sph(mm,  q,  2,  -2,  3,  -2,  3,  N); 

//  zz=zz-Diag  int_A_Four_xy_sph(mm,  q,  2,  0,  1 ,  -2,  3,  N)*Diag_int_A_Four_xy_sph(mm,  q,  2,  -2,  3,  0, 

1 ,  N)/denomin6; 

//b  1  -b2  01 

denomin7=l/Mie2(2,  q,  n_r)+Diag_int_A_Four_xy_sph(mm,  q,  2,  1, 2,  1 , 2,  N); 
zz=zz-Diag_int_A_Four_xy_sph(mm+ 1/2,  q,  2,  0,  1,  1,2,  N)*Diag_int_A_Four_xy_sph(mm-l/2,  q,  2, 

1, 2,  0,  1,  N)/denomin7; 

//bl  -b2  0-1 

denomin8=l/Mie2(2,  q,  n_r)+Diag_int_A_Four_xy_sph(mm,  q,  2,  -1 , 2,  -1 , 2,  N); 
zz=zz-Diag_int_A  Four  xy  sph(mm-l/2,  q,  2,  0,  1,  -1 , 2,  N)*Diag_int_A_Four_xy_sph(mrn+l/2,  q,  2,  - 
1, 2,  0,  1,  N)/denomin8; 

denominh=l/Miel  (1 ,  q+h,  n_r)+Diag_int_A_Four_xy_sph(mm-l ,  q+h,  2,  1 ,  n,  1 ,  n,  N)- 
2*Diag_int_A_Four_xy_sph(mm,  q+h,  2,  -1,  n,  1,  n,  N); 

zz_h=l/Mie2(  1 ,  q^h  ,n_r)  +  Diag_int_A_Four_xy_sph(mm,  q+h,  2,  0,  n,  0,  n,  N)- 
4*OffDiag_int_B_Four_xy_sph(mm-l/2,  q+h,  2,  0,  1,  1 ,  1,  N)A2/denominh; 

denominl  h=l/Miel(2,  q+h,  n_r)+Diag_int_A_Four_xy_sph(mm,  q+h,  2,  0,  2,  0,  2,  N); 
zz_h=zz_h-OffDiag_int_B_Four_xy_sph(mm,  q+h,  2,  0,  1, 0,  2,  N)*OffDiag_int_B_Four_xy_sph(mm, 
q+h,  2,  0,  2,  0,  1,  N)/denominlh; 

//  denomin2h=l/Mie2(3,  q+h,  n_r)+Diag_int_A_Four_xy_sph(mm,  q+h,  2,  0,  3,  0,  3,  N); 

//  zz_h=zz_h-Diag  int_A_Four_xy_sph(mm,  q+h,  2,  0,  1 , 0,  3,  N)*Diag_int_A_Four_xy_sph(mm,  q+h, 

2,  0,  3,  0,  1,  N)/denomin2h; 

denomin3h=l/Miel(2,  q+h,  n  r)+Diag_int_A_Four_xy_sph(mm,  q+h,  2,  2,  2,  2,  2,  N); 
zz_h=zz_h-OffDiag_int_B_Four_xy_sph(mm,  q+h,  2,  0,  1, 2,  2,  N)*OffDiag_int_B_Four_xy  sph(mm, 
q+h,  2,  2,  2,  0,  1 ,  N)/denomin3h; 

denomin4h=l/Miel(2,  q+h,  n_r)+Diag_int_A_Four_xy_sph(mm,  q+h,  2,  -2,  2,  -2,  2,  N); 
zz_h=zz_h-OffDiag_int_B_Four_xy_sph(mm,  q+h,  2,  0,  1,  -2,  2,  N)*OffDiag_int  B  Four  xy  sph(mm, 
q+h,  2,  -2,  2,  0,  1,  N)/denomin4h; 

//  denomin5h=l/Mie2(3,  q+h,  n_r)+Diag_int_A_Four_xy_sph(mm,  q+h,  2,  2,  3,  2,  3,  N); 

//  zz_h=zz_h-Diag_int_A_Four_xy_sph(mm,  q+h,  2,  0,  1 , 2,  3,  N)*Diag_int_A_Four_xy_sph(mm,  q+h, 

2,  2,  3,  0,  1,  N)/denomin5h; 

//  denomin6h=l/Mie2(3,  q+h,  n_r)+Diag_int_A_Four_xy_sph(mm,  q+h,  2,  -2,  3,  -2,  3,  N); 
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//  zz_h=zz_h-Diag_int_A_Four_xy_sph(mm,  q+h,  2,  0,  1,  -2,  3,  N)*Diag_int_A_Four_xy_sph(mm,  q+h, 
2,  -2,  3,  0,  1,  N)/denomin6h; 

//bl-b2  01 

denomm7h=l/Mie2(2,  q+h,  n_r)+Diagint_A_Four_xy_sph(mm,  q+h,  2,  1,2,  1,2,  N); 
zz  h=zz  h-Diag  int  A  Four  xy  sph(mm+l/2,  q+h,  2,  0,  1,  1,2,  N)*Diag  int  A  Four  xy  sph(mm- 
1/2,  q+h,  2,  1 , 2,  0,  1 ,  N)/denomin7h; 
fib  1  -b2  0-1 

denomin8h=l/Mie2(2,  q+h,  n_r)+Diag_int_A_Four_xy_sph(mm,  q+h,  2,  -1, 2,  -1, 2,  N); 
zz_h=zz_h-Diag_int_A_Four_xy_sph(mm- 1/2,  q+h,  2,  0,  1,  -1, 2, 
N)*Diag_int_A_Four_xy_sph(mm+l/2,  q+h,  2,  -1,  2,  0,  1,  N)/denomin8h; 
zz_der  =  (zz_h-zz)/h; 
delt=-zz/zz_der; 
q=q+delt; 
if  (k>20) 
delt=0; 
end 
end 

disp(k); 

disp(-real(q)/(2*imag(q))); 

y=q; 

endfunction 

function  y=VeryStrange(n  r,  N,  q_0) 

//  Solver  for  dipolar  approach  +  off-diagonal  interaction  and  lowest  Mie  resonance  for  closely  packed 
//circular  array  taking  into  account  the  off-diagonal  interaction 
//  n  r  -refractive  index  2.7  -  TiO_2,  3.5  -  GaAs,  2  -  ZnO 
//  n=l  -  dipoles,  m=0  - 1 1 ,  m=  1  - 12,  m=- 1  - 1 
//  N  -  number  of  spheres 

//  mm=0,  1 ,  ...N  angular  momentum  of  mode  mm=N/2  -  most  interesting  case 
//  q_0  -  initial  value  for  iteration  procedure 

//  Assuming  the  distance  between  sphere  is  2,  targets  are  the  quality  factor  and  the  decay  rate  for 

//  a=200nm 

h=0.0000000001; 

q~q_0; 

n=l; 

mm=N/2; 

disp("  1  st  transverse  mode  1 1 "); 

delt  =  1 ; 

k=0; 

while  abs(delt/q)  >  hA(3/2) 
k=k+ 1 ; 

//  bl-al  0-  1-1 

denomin=l/Miel(l,  q,  n_r)+Diag_int_A_Four_xy_sph(mm-l,  q,  2,  1,  n,  1,  n,  N)- 
2*Diag_int_A_Four_xy_sph(mm,  q,  2,  -1,  n,  1,  n,  N); 

zz  =  1/Mie2(l ,  q  ,n_r)  +  Diag_int_A_Four_xy_sph(mm,  q,  2,  0,  n,  0,  n,  N)- 
4*OffDiag_int_B_Four_xy_sph(mm-l/2,  q,  2,  0,  1,  1,  1,  N)A2/denomin  ; 

H bl-a2  00 

denomml  =  l/Miel(2,  q,  n_r)+Diag_int_A_Four_xy_sph(mm,  q,  2,  0,  2,  0,  2,  N); 
zz=zz-OffDiag_int_B_Four_xy_sph(mm,  q,  2,  0,  1, 0,  2,  N)*OffDiag_int_B  Four_xy_sph(mm,  q,  2,  0, 
2,  0,  1 ,  N)/denominl; 

//b  1  -b3  00 

denomm2=l/Mie2(3,  q,  n_r)+Diag_int_A_Four_xy_sph(mm,  q,  2,  0,  3,  0,  3,  N); 
zz=zz-Diag_int_A_Four_xy_sph(mm,  q,  2,  0,  1, 0,  3,  N)*Diag_int_A_Four_xy_sph(mm,  q,  2,  0,  3,  0,  1, 
N)/denomin2; 

Ifb  1  -a2  02 

denomin3=l/Miel(2,  q,  n_r)+Diag_int_A_Four_xy_sph(mm,  q,  2,  2,  2,  2,  2,  N); 
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zz=zz-OffDiag_int  B_Four_xy_sph(mm,  q,  2,  0,  1,  2,  2,  N)*OffDiag_int_B_Four_xy_sph(mm,  q,  2,  2, 

2,  0,  1,  N)/denomin3; 

//bl-a2  0-2 

denomin4=l/Miel(2,  q,  n_r)+Diagjnt_A_Four_xy_sph(mni,  q,  2,  -2,  2,  -2,  2,  N); 
zz=zz-OffDiag_int_B_Four_xy_sph(mm,  q,  2,  0,  1,  -2,  2,  N)*OffDiag_int  B  Four_xy_sph(mm,  q,  2,  -2, 
2,  0,  1,  N)/denomin4; 

//b  1  -b3  02 

denornin5=l/Mie2(3,  q,  n_r)+Diagjnt_A_Four_xy_sph(nrirn,  q,  2,  2,  3,  2,  3,  N); 

zz=zz-Diag_int_A_Four_xy_sph(rnrn,  q,  2,  0,  1,  2,  3,  N)*Diag_int_A_Four_xy_sph(mm,  q,  2,  2,  3,  0,  1, 
N)/denomin5; 

//b  1  -b3  0-2 

denomin6^1/Mie2(3,  q,  n_r)+Diag_int_A_Four_xy_sph(mm,  q,  2,  -2,  3,  -2,  3,  N); 
zz=zz-Diag_int_A_Four_xy  sph(mm,  q,  2,  0,  1,  -2,  3,  N)*Diag_int_A_Four_xy_sph(mm,  q,  2,  -2,  3,  0, 

1 ,  N)/denomin6; 

//b  1  -b2  0 1 

denomin7=l/Mie2(2,  q,  n_r)+Diag_int_A_Four_xy_sph(rrim,  q,  2,  1,2,  1, 2,  N); 

zz=zz-Diag_int  A  Four_xy_sph(rnnri+l/2,  q,  2,  0,  1,  1,  2,  N)*Diag_int_A_Four_xy_sph(mm- 1/2,  q,  2, 

1,  2,  0,  1,  N)/denomin7; 

//b  1  -b2  0-1 

denormn8  1/Mie2(2,  q,  n_r)+Diag_int_A_Four_xy_sph(mm,  q,  2,  -1, 2,  -1, 2,  N); 
zz=zz-Diag_int_A_Four_xy_sph(mm-l/2,  q,  2,  0,  1,  -1, 2,  N)*Diag_int_A_Four_xy_sph(mrn+l/2,  q,  2,  - 
1, 2,  0,  1,  N)/denomin8; 

denominh=l/Miel(l,  q+h,  n  r)+Diag_int_A_Four_xy_sph(rnrn-l,  q+h,  2,  l,n,  l,n,  N)- 
2*Diag_int_A_Four_xy_sph(mm,  q+h,  2,  -1,  n,  1,  n,  N); 

zz_h=l/Mie2(l ,  q+h  ,n_r)  +  Diag _int_A_Four_xy_sph(mm,  q+h,  2,  0,  n,  0,  n,  N)- 
4*OffDiag_int_B_Four_xy_sph(mm-l/2,  q+h,  2,  0,  1,  1,  1,  N)A2/denominh; 

denominlh=l/Miel(2,  q+h,  n_r)+Diag_int_A_Four_xy_sph(mm,  q+h,  2,  0,  2,  0,  2,  N); 
zz  h=zz_h-OffDiag_int_B_Four_xy_sph(mm,  q+h,  2,  0,  1, 0,  2,  N)*OffDiag_int_B  Four_xy_sph(mm, 
q+h,  2,  0,  2,0,  1 ,  N)/denomin  1  h; 

denomin2h  1/Mie2(3,  q+h,  n_r)+Diag_int_A_Four_xy_sph(mm,  q+h,  2,  0,  3,  0,  3,  N); 
zz_h=zz_h-Diag_int_A_Four_xy_sph(mm,  q+h,  2,  0,  1, 0,  3,  N)*Diag_int_A_Four_xy_sph(mm,  q+h,  2, 
0,  3,  0,  1,  N)/denomin2h; 

denomin3h=l/Miel(2,  q+h,  n_r)+Diagjnt_A_Four_xy_sph(mm,  q+h,  2,  2,  2,  2,  2,  N); 
zz_h=zz_h-OffDiagjnt_B_Four_xy_sph(rnm,  q+h,  2,  0,  1, 2,  2,  N)*OffDiag  int  B  Four  xy  sph(mm, 
q+h,  2,  2,  2,  0,  1,  N)/denomin3h, 

denomin4h=l/Miel(2,  q+h,  n_r)+Diag_int_A_Four_xy_sph(mm,  q+h,  2,  -2,  2,  -2,  2,  N); 
zz_h=zz_h-OffDiag_int_B_Four_xy_sph(mm,  q+h,  2,  0,  1,  -2,  2,  N)*OffDiag_int_B_Four_xy_sph(mm, 
q+h,  2,  -2,  2,  0,  1,  N)/denomin4h; 

denorrrin5h=l/Mie2(3,  q+h,  n  r)+Diag  int_A  Four_xy_sph(mm,  q+h,  2,  2,  3,  2,  3,  N); 
zz_h^zz_h-Diag_int_A_Four_xy_sph(mm,  q+h,  2,  0,  1,  2,  3,  N)*Diag_int_A_Four_xy_sph(mm,  q+h,  2, 

2,  3,  0,  1,  N)/denomin5h; 

denomin6h=l/Mie2(3,  q+h,  n_r)+Diag_int_A_Four_xy_sph(mm,  q+h,  2,  -2,  3,  -2,  3,  N); 
zz_h=zz_h-Diag_int_A_Four_xy_sph(mm,  q+h,  2,  0,  1,  -2,  3,  N)*Diag_int_A_Four_xy_sph(mm,  q+h, 

2,  -2,  3,  0,  1,  N)/denomin6h; 

//bl-b2  01 

denomin7h=l/Mie2(2,  q  +  h,  n_r)+Diag_int_A_Four_xy_sph(mm,  q+h,  2,  1,2,  1,2,  N); 
zz_h=zz_h-Diag_int_A_Four_xy_sph(mm+l/2,  q+h,  2,  0,  1,  1,2,  N)*Diag  int_A_Four_xy  sph(mm- 
1/2,  q+h,  2,  1, 2,  0,  1,  N)/denomin7h; 

//b  1  -b2  0-1 

denomin8h=l/Mie2(2,  q+h,  n_r)+Diag_int_A_Four_xy_sph(mm,  q+h,  2,  -1,  2,  -1,  2,  N); 
zz_h=zz_h-Diag_int_A_Four_xy_sph(mm-l/2,  q+h,  2,  0,  1,  -1,  2, 
N)*Diag_int_A_Four_xy_sph(mm+l/2,  q+h,  2,  -1, 2,  0,  1,  N)/denomin8h; 
zz_der  -  ( zz_h-zz)/h; 
delt=-zz/zz_der; 
q=q+delt; 
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//  disp(zz); 
if  (k>20) 
delt=0; 
end 
end 

disp(k); 

disp(-real(q)/(2*imag(q))); 

y=q; 

endfunetion 


function  y=round_besselh(n,  z) 

//  Calculates  Round  Bessel  function  (diverging  wave) 
if  (n~=2)&(n~=l  )&(n~0) 
x  1  =besselh(n+ 1 12,  z)/sqrt(%pi*z); 
else 
if  n==0 

xl  =  -%i*exp(%i*z)/z; 
else 
if  n==l 

xl  =  exp(%i*z)*(-l/z-%i/zA2); 
else 

xl  =  exp(%i*z)*(%i/z-3/zA2-3*%i/zA3); 
end 
end 
end 
y=xl; 

endfunetion 

function  y=round_besselh_sp_der(n,  z) 

//  Calculates  Round  Bessel  function  derivative  (diverging  wave) 
y=-round_besselh(n+l,  z)  +  round_besselh(n,  z)*(n+l)/z; 
endfunetion 

function  y=round_besselj(n,  z) 

//  Calculates  Round  Bessel  function  (standing  wave) 
if  (n“2)&(n“2)&(n~=0) 
x  1  =besselj(n  + 1  /2,  z)/sqrt(%pi*z); 
else 
if  n==0 
x  1  =  sin(z)/z; 
else 
if  n=l 

xl=sin(z)/zA2  -eos(z)/z; 
else 

xl  =  (-sin(z)/z-3*eos(z)/zA2+3*sin(z)/zA3); 
end 
end 
end 
y=xl; 

endfunetion 


function  y=round_besselj_sp_der(n,  z) 

//  Calculates  Round  Bessel  function  derivatives  (standing  wave) 
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y=-round_besselj(n+l,  z)  +  round  bessclj(n,  z)*(n+l)/z; 
endfunction 

function  y=ncatti_besselh(n,  z) 

//  Calculates  Ricatti  Bessel  function  (Hankel) 
x  1  =z*besselh(n+ 1  /2,  z)/sqrt(%pi*z); 
y=xl; 

endfunction 


function  y=ricatti_besselh_der(n,  z) 

//  Calculates  Ricatti  Bessel  function  derivative  (Hankel) 
y=-ricatti_besselh(n+l,  z)  +  ricatti  besselh(n,  z)*(n+l)/z; 
endfunction 

function  y=ricatti_besselj(n,  z) 

//  Calculates  Ricatti  Bessel  function  (Bessel) 
x  1  =z*besselj(n+ 1  /2,  z)/sqrt(%pi*z); 
y=xl; 

endfunction 


function  y=ricatti_besselj_der(n,  z) 

//Calculates  Ricatti  Bessel  function  derivative  (Bessel) 
y=-ricatti_besselj(n+l,  z)  +  ncatti_besselj(n,  z)*(n+l)/z; 
endfunction 


